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Abstract. This paper presents a new approach to verify accuracy of computational simulations. 
We develop mathematical theorems which can serve as robust a posteriori error estimation tech¬ 
niques to identify numerical pollution, check the performance of adaptive meshes, and verify numer¬ 
ical solutions. We demonstrate performance of this methodology on problems from flow thorough 
porous media. However, one can extend it to other models. We construct mathematical properties 
such that the solutions to Darcy and Darcy-Brinkman equations satisfy them. The mathematical 
properties include the total minimum mechanical power, minimum dissipation theorem, reciprocal 
relation, and maximum principle for the vorticity. All the developed theorems have firm mechanical 
bases and are independent of numerical methods. So, these can be utilized for solution verifica¬ 
tion of finite element, finite volume, finite difference, lattice Boltzmann methods and so forth. In 
particular, we show that, for a given set of boundary conditions, Darcy velocity has the minimum 
total mechanical power of all the kinematically admissible vector fields. We also show that a sim¬ 
ilar result holds for Darcy-Brinkman velocity. We then show for a conservative body force, the 
Darcy and Darcy-Brinkman velocities have the minimum total dissipation among their respective 
kinematically admissible vector fields. Using numerical examples, we show that the minimum dissi¬ 
pation and total mechanical power theorems can be utilized to identify pollution errors in numerical 
solutions. The solutions to Darcy and Darcy-Brinkman equations are shown to satisfy a reciprocal 
relation, which has the potential to identify errors in the numerical implementation of boundary 
conditions. It is also shown that the vorticity under both steady and transient Darcy-Brinkman 
equations satisfy maximum principles if the body force is conservative and the permeability is ho¬ 
mogeneous and isotropic. A discussion on the nature of vorticity under steady and transient Darcy 
equations is also presented. Using several numerical examples, we will demonstrate the predictive 
capabilities of the proposed a posteriori techniques in assessing the accuracy of numerical solutions 
for a general class of problems, which could involve complex domains and general computational 
grids. 


1. INTRODUCTION 

1.1. Validation and verification (V&V). Errors can arise in both physical modeling and 
numerical simulation. The study of errors due to physical modeling is referred to as validation, and 
the study of error in a numerical simulation is referred to as verification. As Blottner [Blottner, 
1990] nicely puts it, validation is to solve right governing equations'''’ and verification is to solve 
governing equation right”. Validation errors arise when a model is used out of its application 

Key words and phrases, accuracy assessment; validation and verification (V&V); Darcy model; Darcy-Brinkman 
model; mechanical dissipation; reciprocal relations; vorticity; maximum principles. 
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range. The errors in the verification, on the other hand, can arise from three broad sources includ¬ 
ing numerical errors, round-off errors (due to the finite precision arithmetic), and programming 
mistakes [Oberkampf et ah, 2004]. Basically, the verification is to ensure that the code produces a 
solution with some degree of accuracy, and the numerical solution is consistent. Verification itself is 
conducted into two modes; verification of code and verification of calculation [Roache, 1998; Rider 
et ah, 2016]. Verification of code addresses the question of whether the numerical algorithms have 
been programmed and implemented correctly in the code. The two currently popular approaches 
to verify a code are the method of exact solutions (MES) and the method of manufactured solutions 
(MMS). More thorough discussions on MES and MMS can be found in [Knupp and Salari, 2003; 
Roy et ah, 2004; Roache, 1998]. 

Verification of calculation (which is also referred to as solution verification) estimates the overall 
magnitude (not just the order) of the numerical errors in a calculation, and the procedure invariably 
involves a posteriori error estimation [Salari and Knupp, 2000]. The numerical errors in the solution 
verihcation can arise from two different sources including discretization errors and solution errors. 
The discretization errors refer to all the errors caused by conversion of the governing equations 
(PDEs and boundary conditions) into discrete algebraic equations whereas the solution errors refer 
to the errors in approximate solution of the discrete equations. The numerical errors may arise from 
insufficient mesh resolution, improper selection of time-step, and incomplete iterative convergence. 
Eor more details on verification of calculation, see [Roy et ah, 2004; Salari and Knupp, 2000; Roache, 
1997, 1998; Babuska and Oden, 2004; Oberkampf and Blottner, 1998; Oberkampf et ah, 2004; Rider 
et ah, 2016]. 

1.2. A posteriori techniques. The aim of a posteriori error estimation is to assess the 
accuracy of the numerical approximation in the terms of known quantities such as geometrical 
properties of computational grid, the input data, and the numerical solution. A posteriori error 
techniques monitor various forms of the error in the numerical solution such as velocity, stress, 
mean fluxes, and drag and lift coefficients [Becker and Rannacher, 2001]. Such error estimation 
differ from a priori error estimates in that the error controlling parameters depend on unknown 
quantities. A priori error estimation investigates the stability and convergence of a solver and can 
give rough information on the asymptotic behavior of errors in calculations when grid parameters 
are changed appropriately [Ainsworth and Oden, 1997]. 

Moreover, discontinuities and singularities commonly occur in fluid dynamics, solid mechanics, 
and structural dynamics and it is a subject of intense research and grave concern in error estimation 
[Oberkampf et ah, 2004]. Some representative works on pollution errors due to singularities and 
discontinuities are [Babuska and Oh, 1987; Babuska et ah, 1995, 1997; Oden et ah, 1998; Roache, 
1998; Botella and Peyret, 2001]. 

The question that still remains is how to verify the accuracy of numerical solutions, especially 
for realistic problems. Below are some specific challenges that a computational scientist may face 
in using numerical simulators: 

(i) How much mesh refinement is required to obtain solutions with desired degree of accuracy for 
a problem that does not have an analytical or reference solution? 

(ii) Will the chosen mesh be able to resolve singularities in the solution and avoid pollution errors? 
That is, can we identify whether a particular type of mesh suffers from pollution errors for a 
problem with singularities? 

(iii) Has the computer implementation been done properly? 
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(iv) Is the chosen numerical formulation accurate/appropriate for the chosen problem? 

In the literature, one finds the usual approach of employing tools from functional analysis to obtain 
a priori estimates and assess stability. For example, see [Brezzi and Fortin, 1991; Babuska and 
Strouboulis, 2001]. On the other hand, this paper aims to address the aforementioned challenges 
by providing various a posteriori techniques with firm mechanics underpinning. 

In the current study, we develop new methodology to verify the accuracy and convergence 
of numerical solutions. In addition, the presented a posteriori error estimation approach is able 
to identify pollution in computational domains and check whether adaptive grids can resolve the 
numerical pollution. The proposed criteria is independent of employed numerical method. It means, 
our technique is able to verify the solution of all numerical methods including finite element, hnite 
volume, finite difference and lattice Boltzmann. Hence, we named it mechanics-based solution 
verihcation. Then, we show performance of this powerful tool on problems from flow thorough 
porous solid. But, it is not restricted to porous media models and it is extendable to other problems. 


1.3. Popular porous media models. The simplest and yet the most popular model that 
describes the flow of an incompressible fluid in a rigid porous media is the Darcy model [Darcy, 
1856]. The Darcy equations posed several challenges to the computational community but played a 
crucial role in the development of mixed and stabilized formulations in finite element method (FEM) 
[Masud and Hughes, 2002; Nakshatrala et ah, 2006]. Brinkman [Brinkman, 1947a,b] proposed an 
extension to the Darcy model which is commonly referred to as the Darcy-Brinkman model. In 
addition to drag between the fluid and porous solid, Darcy-Brinkman model accounts for viscous 
term (friction between fluid layers). It is, however, not possible to obtain analytical solutions under 
these models, and one commonly seeks numerical solutions for realistic problems. 

Since the aim of this paper is a posteriori error estimation, we assume that the code has already 
been verihed for the Darcy and Darcy-Brinkman class of problems so that programming mistakes 
are not an issue. Likewise, we are not also concerned with validation. It means that the Darcy and 
Darcy-Brinkman models are physically adequate to model the problems. Moreover, the solution 
verihcation requires conhrmation of grid convergence which is one of the most common and reliable 
error estimation methods [Roache, 1997]. Similar to grid convergence studies, we address only 
solution verihcation. 


1.4. Organization of the paper. Section 2 presents the governing equations arising from 
the Darcy and Darcy-Brinkman models. In Section 3, we propose various mathematical properties 
that the solutions to these governing equations satisfy. We also discuss how these properties can 
be utilized as robust a posteriori criteria to assess the accuracy of numerical solutions. Section 
4 presents several steady-state numerical results to illustrate the predictive capabilities of the 
proposed a posteriori criteria with respect to singularities, pollution errors, and discretization 
errors in the implementation of (Neumann) boundary conditions. In Section 5, we utilize synthetic 
reservoir data to demonstrate the usefulness of the proposed a posteriori techniques, especially for 
problems involving spatially heterogeneous permeability properties. Section 6 discusses a posteriori 
criteria for transient problems, and presents representative numerical results in support of the 
theoretical predictions. Finally, conclusions are drawn in Section 7. 

3 



2. DARCY AND DARCY-BRINKMAN MODELS 


Let ri C be an open and bounded domain, where “nd” denotes the number of spatial 
dimensions. We shall denote the set closure of Q by fl. Let <9^ := denote the boundary, 

which is assumed to be piecewise smooth. A spatial point in Q is denoted by x. The spatial gradient 
and divergence operators are, respectively, denoted as grad[-] and div[-]. Let v : —)■ denote 

the velocity field and p : —)• M denote the pressure field. The symmetric part of the gradient of 
velocity is denoted by D(x). That is, 

D(x) = ^ (grad[v] + grad[v]'^) (2.1) 

The unit outward normal to the boundary is denoted as n(x). The boundary is divided into two 
parts: T’' and Th T’^ is the part of the boundary on which the velocity is prescribed, and T* is 
that part of the boundary on which the traction is prescribed. For mathematical well-posedness, 
we have T’' n T* = 0 and T’^ U T* = 

The porous media models that will be considered in this paper are the Darcy and Darcy- 
Brinkman models. Both these models describe the flow of an incompressible fluid through rigid 
porous media. We completely neglect the motion of the porous solid. The Cauchy stress in the 
Darcy and Darcy-Brinkman models, respectively, take the following form: 

T(x) = —p(x)I (2.2a) 

T(x) = —p(x)I + 2|uD(x) (2.2b) 

where I denotes the second-order identity tensor, and g, is the dynamic coefficient of viscosity. The 
steady-state governing equations based on the Darcy model can be written as follows: 


q;(x)v(x) -|- grad[p(x)] = /9b(x) 

in D 

(2.3a) 

div[v(x)] = 0 

in D 

(2.3b) 

v(x) • n(x) = Un(x) 

on 

(2.3c) 

p(x) = po(x) 

on T* 

(2.3d) 


where a(x) is the drag coefficient, p is the density of the fluid, b(x) is the specific body force, 
Un(x) is the prescribed normal component of the velocity, and po(x) is the prescribed pressure. The 
steady-state governing equations based on the Darcy-Brinkman model take the following form: 


a(x)v(x) -|- grad[p(x)] — div [2/iD] = pb(x) 

in D 

(2.4a) 

div[v(x)] = 0 

in D 

(2.4b) 

v(x) = vP(x) 

on r 

(2.4c) 

TS(x) = tP(x) 

on r* 

(2.4d) 


where vP(x) is the prescribed velocity vector, and tP(x) is the prescribed traction. We shall call a 
vector-held to be Darcy velocity if satishes equations (2.3a)-(2.3d). We shall call a vector held to 
be Darcy-Brinkman velocity if it satishes equations (2.4a)“(2.4d). 

Equations (2.3a) and (2.4a) can be obtained from the balance of linear momentum under the 
mathematical framework offered by the theory of interacting continua [Nakshatrala and Rajagopal, 
2011]. The drag term q;(x)v(x) models the frictional force between the huid and the porous 
solid. The term div[2^D] in the Darcy-Brinkman model arises due to the internal friction between 
the layers of the huid. The pressure p(x) is an undetermined multiplier that arises due to the 
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enforcement of the incompressibility constraint given by equations (2.3b) and (2.4b). The drag 
coefficient is related to the coefficient of viscosity of the fluid and the permeability k{x.) as follows: 


a(x) 


/c(x) 


(2.5) 


In general, it is not possible to obtain analytical solutions to the systems of equations given 
by either (2.3a)-(2.3d) or (2.4a)-(2.4d). Hence, one needs to resort to numerical solutions. This 
paper does not coneern with developing new numerical formulations to solve the aforementioned 
mathematical models. The paper instead focuses on deriving mathematical properties with firm 
mechanics underpinning that the solutions to these mathematical models satisfy. We shall then 
illustrate how these mathematical properties can serve as robust “a posteriori” criteria to assess the 
accuracy of numerical solutions. 


3. MATHEMATICAL PROPERTIES: STATEMENTS AND DERIVATIONS 

In the remainder of this paper, we shall refer to a vector field v : H —)• as kinematically 

admissible if it satisfies the following conditions; 

(i) v(x) is solenoidal (i.e., div[v(x)] = 0 in H), and 

(ii) v(x) satisfies the boundary conditions. 

It needs to emphasized that a kinematically admissible vector field need not satisfy the balance of 
linear momentum given by equation (2.3a) or (2.4a). Clearly, the Darcy velocity and the Darcy- 
Brinkman velocity are kinematically admissible vector fields. For some of the results presented 
in this paper, we will need the body force to be conservative, which is a common terminology in 
potential theory [Kellogg, 2010]. The body force /9b(x) is said to be conservative if there exists a 
scalar potential V’(x) such that pb(x) = —grad[V']. We shall define the dissipation functional as 
follows: 


J a(x)v(x) • v(x) dD Darcy model 

\ a(x)v(x) ■ v(x) dD + 2/uD(x) ■ D(x) dD Darcy-Brinkman model 


Since a > 0 and ^ > 0, it is straightforward to check that 4>[v] is a norm. In fact, it can be shown 
that <I>[v] under the Darcy model is equivalent to the natural norm in (L^(D))"^'^, where (L^(D))"''^ 
is a space of square integrable vector fields defined from D to Similarly, it can be shown that 
4>[v] under the Darcy-Brinkman model is equivalent to the natural norm in which is a 

Sobolev space. For further details on function spaces and norms, refer to [Evans, 1998]. 

In this section, we shall present four important mathematical properties that the solutions 
to Darcy equations and Darcy-Brinkman equations satisfy. These properties will be referred to 
as (i) the minimum total mechanical power theorem, (ii) the minimum dissipation theorem, (hi) 
reciprocal relation, and (iv) maximum principle for vorticity. As a passing comment, we will 
employ the minimum dissipation theorem to show the uniqueness of solution for Darcy equations 
and Darcy-Brinkman equations. In the porous media literature, these results have neither been 
discussed nor utilized to solve problems. More importantly, these results have not been used to 
assess the accuracy and convergence of numerical solutions of porous media models. For example, 
it will be shown that the minimum total mechanical power theorem can be utilized to assess the 
accuracy of the implementation of both Dirichlet and Neumann boundary conditions. On the other 
hand, the minimum dissipation theorem can be utilized to identify pollution errors in the numerical 
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solution. Herein, we give detailed mathematical proofs for Darcy-Brinkman equations. We however 
provide comments on the corresponding proofs for Darcy equations. 

Theorem 1 (Minimum total mechanical power theorem). Let v(x) be the Darcy-Brinkman ve¬ 
locity vector field. Then, any kinematically admissible vector field v(x) satisfies the following in¬ 
equality: 


eTMp[v] < £tmp[v] 


(3.2) 


where 

eTMp[z] := - [ /C'b(x) • z(x) dO. - [ tP(x) • z(x) dD (3.3) 

^ Jn Jr* 

That is, for given boundary conditions, body force and tractions; the Darcy-Brinkman velocity will 
have the minimum total mechanical power among all the possible kinematically admissible vector 
fields. 

Proof. Let 


5v(x) := v(x) — v(x) 

(5D(x) := D(x) — D(x) 

From the hypothesis of the theorem, (5v(x) satisfies the following relations: 

(5v(x) = 0 Vx G dLl 
div[(5v] = 0 Vx G D 

Let us start with the dissipation due to the vector field v(x): 


(3.4a) 

(3.4b) 


(3.5a) 

(3.5b) 


‘h[v(x)] := / a(x)v(x) • v(x) dD + / 2^D(x) • D(x) dD 

Jn Jn 

= f q;(x) ((5v(x) + v(x)) • ((5v(x) + v(x)) dD + f 2^ (5D(x) + D(x)) • ((5D(x) + Dx)) dD 

Jn Jn 

>2 / a(x)(5v(x) • v(x) dD + 2 f 2^5D(x) • D(x) dD + 4>[v(x)] (3.6) 

Jn Jn 

Using equations (2.4a) and (2.2b) the first integral in the above equation can be written as follows: 

f q;(x)5v(x) • v(x) dD = f 5v(x) • (pb(x) + div [T(x)]) dD (3.7) 

J ^ J ^ 

The symmetry of D(x) allows the second integral to be written as follows: 

f 2fi5D{x) ■ D(x) dD = f 2p. grad[(5v(x)] • D(x) dD 

Jn Jn 

= f grad[(5v(x)] • (T(x) +p(x)I) dD 

Jn 

= f grad[(5v(x)] • T(x) dD + / div[(5v(x)] • p(x) dD 

Jn Jn 
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Noting that div[(5v] = 0 in O and by employing Green’s identity, we have the following: 

f 2//(5D(x) • D(x) dn = / div[T^(x) (5v(x)] dfi — / 5v(x) • div[T(x)] dO 

Jo, Jq Jq, 

= [ 5v{x) ■ tP(x) dr — / (5v(x) • div[T(x)] dfl 

ir* Jn 

From equations (3.7) and (3.8), inequality (3.6) can be written as follows; 

^[v(x)] > $[v(x)] + 2 / (5v(x) • /3b(x) dfl + 2 f 6v{x) ■ tP(x) dr 

Jn Jr* 

This completes the proof. 


(3.8) 

(3.9) 
□ 


Remark 1. For Darcy equations (2.3a)-(2.3d), one can state the minimum total mechanical 
power theorem as follows: For given boundary conditions, and body force; the Darcy velocity, v(x), 
has the minimum total mechanical power among all the kinematically admissible vector fields. That 
is, 

- [ v(x) • /ob(x) dD+ [ po(x)v(x) • n(x) dF < ;^<h[v] - [ v(x) • pb(x) dR 

2 Jn Jr* 2 Jn 

+ / po(x)v(x) • n(x) dr Vv(x) (3.10) 
Jr* 

Theorem 2 (Minimum dissipation inequality). Let v(x) be the Darcy-Brinkman velocity vector 
field, and let F’’ = dLl. Then, any kinematically admissible vector field v(x) satisfies the following 
inequality: 


<I>[v] < ^[v] (3-11) 

That is, for given velocity boundary conditions and conservative body force, the Darcy-Brinkman 
velocity has the minimum total dissipation due to drag and internal friction of all the possible 
kinematically admissible vector fields. 

Proof. We shall employ the notation introduced in equation (3.4). Recall that the mechanical 
dissipation under the Darcy-Brinkman model is 

= [ C(dv{x) ■ v(x) dD -|- f 2//(5D(x) • D(x) dR (3-12) 

Jn Jn 

Let us start with the difference in total dissipation, and from inequality (3.6) we have: 


<l>[v(x)] — <l>[v(x)] > 2 / a(5v(x) • v(x) dD + 2 / 2|U(5D(x) • D(x) dR (3.13) 

Jn Jn 

Using Green’s identity, the first integral can be simplified as follows: 

f a5v(x) • v(x) dD = f (5v(x) • grad['0(x) — p(x)] dD-|- f (5v(x) • div[2/rD(x)] dD 

= f (5v(x) • n(x) (V’(x) — p{x)) dr — / div[(5v(x)] ('0(x) — p{x)) dD 

Jdn Jn 

J- f (5v(x) • div[2|uD(x)] dD (3-14) 

Jn 
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Noting the symmetry of D(x) and using Green’s identity, the total dissipation due to internal 
friction can be simplified as follows: 


f 2^5D(x) • D(x) dn = / 2|U grad[(5v(x)] • D(x) dO 

J Q. Jft 

= f 2/r(5v(x) • D(x)n(x) dT — f (5v(x) • div[2^D(x)] dfi (3.15) 

Jdn Jn 

From equations (3.7)-(3.8), the total dissipation due to drag and internal friction satisfies; 

<l>[v(x)] — <l>[v(x)] > 2 f (5v(x) • n(x) ('0(x) — p(x)) dF — 2 / div[(5v(x)] (^/)(x) — p(x)) dfl 
Jan Jn 

+ 2 / 2;U(5v(x) • D(x)n(x) dF = 0 (3.16) 

Jan 

This completes the proof. □ 


The solutions to Darcy and Darcy-Brinkman equations posses reciprocal relations similar to the 
famous Betti’s reciprocal relation in the theory of linear elasticity [Truesdell and Noll, 2004; Sadd, 
2009] and to a classical reciprocal relation in the area of creeping flows [Guazzelli and Morris, 
2012]. The Betti’s reciprocal relation is often employed to solve a class of problems in linear 
elasticity, which otherwise may be difficult to solve. Mathematically, the Betti’s reciprocal relation 
is equivalent to the existence and symmetry of Green’s function. We now precisely state a reciprocal 
relation that the solutions of Darcy-Brinkman equations satisfy, and then provide a mathematical 
proof. 


Theorem 3 (Reciprocal relation). Assume that vP(x) = 0 on F^. Let {vi(x),pi(x)} and 
{v 2 (x),p 2 (x)} be the solutions of equations (2.4a)-(2.4d) for the prescribed data {bi(x),tj’(x)} 
and {b 2 (x), t 2 (x)}, respectively. Then, these fields satisfy the following relation: 

/ yobi(x) • V 2 (x) dD -I- f (x) • V 2 (x) dF = f ph2{x.) ■ vfix) dLl + [ t^(x) • vi(x) dF (3.17) 
Jn Jn Jr* 

Proof. Let us start with the left side of equation (3.17). Noting that vP(x) = 0 on F^ and 
F’’ U F* = dn, one can proceed as follows: 

[ /9bi(x) • V 2 (x) dLl+ [ t^(x) • V 2 (x) dr = / /9bi(x) • V 2 (x) dLl+ [ (Tin(x)) • V 2 (x) dF 

Jn Jr* Jn Jr* 

= [ /9bi(x) • V 2 (x) dLl+ [ (Tin(x)) • V 2 (x) dF 

Jn Jan 

= [ /9bi(x) • V 2 (x) dLl+ [ div [T^(x)v 2 (x)] dD 

Jn Jn 

= f (/9bi(x) -I- div[Ti]) • V 2 (x) dD -I- / Ti(x) • grad [V 2 ] dD 
Jn Jn 

= f q;(x)vi(x) • V 2 (x) dLl+ [ (-pi(x)I -|- 2/iDi(x)) • grad [V 2 ] dD 
Jn J n 

= / q;(x)vi(x) • V2(x) dD — / pi(x)div[v2] dD-|- / 2^Di(x) • D2(x) dD 

Jq Jq Jfl 



In the above step, we have used the fact that Di(x) is a symmetric second-order tensor. Since 
div [V 2 ] = 0 we have the following: 


[ / 9 bi(x) • V2(x) d^l+ [ tj’(x) • V2(x) dr = / a(x)vi(x) • V2(x) dil. + [ 2 ^Di(x) • D2(x) dll 

Jn Jr^ Jn Jn 

Similarly, it can be shown that the right side of equation (3.17) is also equal to 

f a(x)vi(x) • V2(x) dQ+ [ 2 ^Di(x) • D2(x) dO 

Jfi Jft 

This completes the proof. 


□ 


The following notation will be used later to verify the reciprocal relation: 

left integral — right integral 

^reciprocal •— \ 7 ~. 7 \ 

left integral 

where the left and right integrals are, respectively, defined as follows: 


left integral := 
right integral := 


/ 

Jn 

I 

Jn 


/jbi(x) • V 2 (x) dn + 
/9b2(x) • vi(x) dn + 



ti(x) • V 2 (x) dr 
t2(x) • Vi(x) dr 


(3.18) 


(3.19a) 

(3.19b) 


Remark 2. The corresponding reeiproeal relation for Darey equations can he written as follows: 
Assume that Un(x) = 0 on r’'. Let {vi(x),pi(x)} and {v 2 (x),p 2 (x)} he the solutions of equations 
(2.4a)-(2.4d) for the prescribed data {bi(x),poi(x)} and {b 2 (x),po 2 (^)}; respeetively. Then, these 
fields satisfy the following relation: 

[ pbi(x) • V2(x) dll - [ poi(x)n(x) • V2(x) dr 

Jn Jr* 

= [ pb 2 (x) • Vi(x) dLl- [ P 02 (x)n(x) • vi(x) dV (3.20) 

Jn Jr* 

Remark 3. It needs to he emphasized that the reeiproeal relation will not directly he able to 
assess the aecuraey of the pressure field in the computational domain. The reeiproeal relation is 
ideal for assessing the aceuraey of the velocity vector field in the domain and the aceuraey of the 
implementation of prescribed traction boundary conditions. This reciprocal relation will not be able 
to provide information about the aeeuraey of the implementation of non-zero velocity boundary 
conditions. 


Next, we discuss in the form of a theorem on the nature of the vorticity under the Darcy and 
Darcy-Brinkman models. To this end, 

a;(x) = curl[v(x)] (3-21) 

In a Cartesian coordinate system, the components of vorticity take the following form: 

_ ^ 

dy dz ’ dz dx ’ dx dy 

It should be emphasized that curl[-] operator is defined only in (i-e., the three-dimensional 
Euclidean space). However, for two-dimensional problems, one can consider the vorticity as follows: 

cj(x, y) = W;j(x, y)e^ (3.23) 

where z denotes the axis perpendicular to the two-dimensional plane in which the problem is 
defined, and e^ is the unit vector along the z-direction. 
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Theorem 4 (On the nature of vorticity under Darcy and Darcy-Brinkman models). Assume that 
the medium is isotropic and homogeneous (i.e., a(x) is a constant scalar), the body force is a 
conservative vector field (i.e., /ob(x) = —grad[r/)(x)] and the response is steady-state. Then the 
vorticity vanishes under Darcy equations. Under Darcy-Brinkman equations, the vorticity is an 
eigenvector of the Laplacian with 1/k as the eigenvalue. Moreover, the vorticity satisfies a maximum 
principle in which the non-negative maximum and the non-positive minimum occur on the boundary. 


Proof. By taking the curl on both sides of the balance of linear momentum under Darcy 
equations (2.3a) we obtain the following: 

curl[av] = —curl[grad[V' + p]] = 0 (3.24) 

Since a is spatially homogeneous scalar, one can conclude that the vorticity vanishes under the 
Darcy model. 

The incompressibility constraint implies that the balance of linear momentum under the Darcy- 
Brinkman model can be written as follows: 


q;v(x) -|- grad[p(x)] — ;udiv[grad[v]] = —grad[V’] (3.25) 

By taking curl on both sides of the above equation, we get the following: 

Aa;(x) = -cu(x) (3.26) 

k 

where A denotes the Laplacian operator and k is the permeability. The above equation is a vector 
eigenvalue problem in which the vorticity vector is the eigenvector, and 1/A: is the corresponding 
eigenvalue. For two-dimensional problems, we have the following: 

A2DW2(x,y) = ia;2(x,y) (3.27) 

where A 2 D denotes the two-dimensional Laplacian operator. The above equation is a scalar eigen¬ 
value problem in which w^(x, y) is the eigenvector and 1/A is the corresponding eigenvalue. 

Since 1/k > 0, this equation is commonly referred to as diffusion with decay, which is a linear 
self-adjoint elliptic partial differential equation. It is well-known that such a partial differential 
equation satisfies a maximum principle [Gilbarg and Trudinger, 2001]. Mathematically, the max¬ 
imum principle for the vorticity under the Darcy-Brinkman model takes the following form: If 
Wzi'x.) G C‘^{D) n (7^(11), then the non-negative maximum and the non-positive minimum occur on 
the boundary. That is, 


max[a;^(x)] < max 0, maxa;^(x) 
xen L xe90 


(3.28a) 


min[ti; 2 (x)] > min 0, min cj 2 (x) 
xeo L xean 


(3.28b) 


where C‘^{Q) denotes the set of twice differentiable functions defined on fl, and C''^(n) is the set of 
functions that are continuous to the boundary. □ 


The maximum principle is used to obtain physical meaningful numerical solutions (see [Naksha- 
trala et ah, 2015; Mudunuru and Nakshatrala, 2016] and References therein). Herein, it can be 
utilized to verify the accuracy of numerical solutions by plotting vorticity and checking whether the 
non-negative maximum and non-positive minimum of the vorticity occur on the boundary. For the 
Darcy model with isotropic and homogeneous medium properties and conservative body force, it 
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can be shown that the vorticity vanishes (i.e., <^(x) = 0). However, it should be noted that hetero¬ 
geneity, pressure-dependent viscosity, or non-conservative body force can introduce vorticity under 
the Darcy model. All the above results can serve as invaluable tools to assess the performance of a 
numerical formulation to verify a computer implementation, and to provide metrics for numerical 
convergence. 

4. STEADY-STATE NUMERICAL RESULTS 

We shall first non-dimensionalize the governing equations by choosing primary variables that 
seem appropriate for problems arising in modeling of flows through porous media. This non- 
dimensional procedure is different from the standard non-dimensionalization procedure for incom¬ 
pressible Navier-Stokes in the choice of primary variables. In the standard non-dimensionalization 
of Navier-Stokes equations, one employs characteristic velocity v, characteristic length L and den¬ 
sity of the fluid p as primary variables. We shall choose L (reference length in the problem), g 
(acceleration due to gravity) and Patm (atmospheric pressure) as the reference quantities. Using 
these reference quantities, we define the following non-dimensional quantities; 

_x_ v„ vP— T — b_ p 

X = — V = —vP = —T = - b = — V — - 

^ V 9^ V 9^ Patm 9 Patm 

Po - P9L _ _ P^fgjL - 

Po = -, P = -, a =-, P = -, k = 

Patm Patm Patm Patm 

where the non-dimensional quantities are denoted by superposed bars. The gradient and divergence 
operators with respect to x are denoted as grad[-] and div[-], respectively. The scaled domain Dgcaied 
is defined as follows: A point in space with position vector x G Dgcaied corresponds to the same 
point with position vector given by x = xL G D. Similarly, one can define the scaled boundaries: 
SDscaled, and Using the above non-dimensionalization procedure, Darcy-Brinkman 


equations can be written as follows: 

a(x)v(x) -h grad[p(x)] - div[2pD(x)] = pb(x) in ^scaled (4.2a) 

div[v(x)] =0 in ^scaled (4-2b) 

v(x)=vP(x) onU^glgd (4.2c) 

t(x)=F(x) onr*^gigd (4.2d) 


Similarly, one could write the corresponding non-dimensional form of Darcy equations. For sim¬ 
plicity, the “over-lines” will be dropped in the remainder of the paper. 

We shall use several flow through porous media problems with different boundary conditions to 
illustrate that the proposed a posteriori techniques can be used as good measures for the accuracy 
and convergence of numerical results. We have employed P2P1 (which is based on six-node triangle 
element T6) and Q2P1 (which is based on nine-node quadrilateral element Q9) interpolations 
available in COMSOL [COM, 2013]. Consistent SUPG stabilization has been employed if the 
interpolations (i.e., Q2P1 and P2P1) violate the LBB inf-sup stability condition [Hughes et ah, 
1986; COM, 2013]. Typical structured finite elements utilized in this paper are shown in Figure 
1. Unless mentioned otherwise, all the elements in a quadrilateral mesh are squares and all the 
elements in a triangular mesh are right-angled isosceles triangles. In this numerical solution study, 
we have taken h (the maximum element size) to be equal to the length of the side for square 
elements, and to the length of the base (or height) for right-angled isosceles triangles. 
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Table 1. Body force problem: Non-dimensional parameters used in the problem. 


Parameter 

Value 

a 

1 

T 

1 and 0.001 

P 

1 

L 

1 

b(x) 

a[sin(7rx), cos(7ry)] 


4.1. Body force problem. The test problem is pictorially described in Figure 2. The non- 
dimensional parameters used in the numerical simulation are provided in Table 1. The conservative 
body force is taken as /ob(x) = 10[sin(7rx), cos(7ry)] (i.e., a = 10). Velocity boundary condition is 
prescribed on the entire boundary (i.e., F^ = dO,). For the Darcy-Brinkman model, we assume 
vP(x) = 0; and for the Darcy model, we assume Vn(x) = 0. Figure 3 shows the minimum total 
mechanical power and the minimum dissipation with mesh refinement. The numerical results for 
the reciprocal relation with mesh refinement are shown in Figure 4. All the numerical results are 
in accordance with the theoretical predictions for this test problem. 

4.2. Lid-driven cavity problem. The two-dimensional lid-driven cavity problem is a bench¬ 
mark study widely used to investigate the accuracy of numerical formulations for various fluid mod¬ 
els [Burggraf, 1966; Erturk, 2009]. Figure 5 provides a pictorial description of the problem. The 
domain of the problem is a bi-unit square. Velocity boundary condition is prescribed on the entire 
boundary (i.e., F^ = dD) which implies that the minimum dissipation theorem is also applicable 
to this problem. It should be noted that lid-driven cavity problem is not compatible with Darcy 
equations which need only the normal component of the velocity to be prescribed on the boundary. 
However, the lid-driven cavity problem demands the prescription of the entire velocity vector field 
on the boundary. We shall therefore employ Darcy-Brinkman model in our numerical simulations. 

It is crucial to note that the solution to the lid-driven cavity problem has singularities at the 
top corners which arise due to velocity discontinuity on the boundary [Botella and Peyret, 1998; 
Batchelor, 2000]. 

Herein, we use an adaptive mesh to resolve the singularities in the solution. The non-dimensional 
parameters used in the lid-driven cavity problem are presented in Table 2. Figures 6(a) and 
6(b) show the uniform structured mesh and the adaptive mesh, respectively. For the adaptive 
mesh, we generate fine grid in the region close to top lid (see Figure 6(b)). Figures 6(c) and 6(d) 
show variation of the minimum dissipation with mesh refinement for respective grids. Under the 
adaptive mesh, the dissipation decreased uniformly and converged to a constant value with mesh 
refinement. On the other hand, the dissipation increased monotonically with mesh refinement under 
the structured mesh. More importantly, due to the presence of singularities and pollution error, 
the dissipation did not hit a plateau even for very fine meshes (i.e., Ijh > 220). The dissipation 
reached a plateau relatively quickly under the adaptive mesh (say lfh = 10). This problem clearly 
illustrates that the minimum dissipation theorem can he used to identify pollution errors and to 
check performance of adaptive mesh in numerical solutions which is one of the main findings of 
this paper. 
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Table 2. Lid-driven cavity problem: Non-dimensional parameters used in the problem. 


Parameter 

Value 

a 

I 

V 

I 

P 

I 

L 

I 


4.3. Pipe bend problem. As another application of the proposed techniques, we consider 
an engineering problem commonly found in the fluid mechanics literature, which is called the pipe 
bend problem (for example, see References [Borrvall and Petersson, 2003; Hansen et ah, 2005; Aage 
et ah, 2008; Pingen et ah, 2009; Challis and Guest, 2009; Hassine, 2012]). In the pipe bend problem, 
the computational domain is a square with L = 1. 

4.3.1. Velocity boundary condition. The problem is pictorially described in Figure 7. An inflow 

parabolic velocity is enforced on a part of the left boundary and an outflow parabolic velocity 
on a segment of the bottom. Each parabolic velocity profile has a unit maximum value (i.e., 
'^Xmax = 1 or = 1). Elsewhere, the homogeneous velocity is prescribed (i.e., vP(x) = 0 for 

the Darcy-Brinkman model and Un(x) = 0 for the Darcy model). The velocity boundary condition 
makes the problem compatible with the total mechanical power and the dissipation theorems. The 
non-dimensional parameters used in the problem are presented in Table 3. Eigure 8 shows the 
variation of minimum total mechanical power and minimum dissipation with mesh refinement for 
the quadrilateral and triangular elements. The result of the numerical solutions verification are 
presented in Eigure 8(a) and 8(b) for the minimum dissipation and the total mechanical power, 
respectively. The numerical error decreases and converges uniformly. 

4.3.2. Parabolic velocity-pressure boundary condition. A pictorial description of the problem is 

given by Eigure 9. The traction is prescribed on a part of the bottom boundary (i.e., tP(x) = 
—Patmn(x) on r*). The velocity has parabolic profile with unit maximum value (i.e., = 1) 

prescribed on a segment of the left boundary. Elsewhere, the homogeneous velocity is enforced. 
On account of the traction and parabolic velocity boundary conditions, current problem is not 
compatible with the minimum dissipation and reciprocal theorems. It is only compatible with 
the total mechanical power theorem. It should be noted that the solution to the problem has 
the singularity at near corners of the outlet (i.e., E*). Hence, we again use an adaptive mesh to 
resolve the pollution in the solution. The non-dimensional parameters using in the problem are 
provided in Table 3. Eigure 10 depicts the variation of the minimum total mechanical power with 
mesh refinement for the quadrilateral and triangular elements. The top figures show the uniform 
structured and adaptive meshes. Under the structured mesh, the total mechanical power increased 
uniformly with mesh refinement due to the polluted area in the computational domain. The error 
decreased uniformly and converged to a constant value with mesh refinement using the adaptive 
mesh. So, in addition to error estimation capability, the total mechanical power theorem can 
be used to identify the pollution errors in the numerical solutions and check the performance of 
adaptive mesh. 

4.4. Pressure slab problem. Eigure II provides a pictorial description of the problem. The 
non-dimensional parameters used in the numerical simulation are provided in Table 4. The domain 
is a kE X L rectangle. The homogeneous velocity boundary condition is enforced on the top and 
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Table 3. Pipe bend problem: Non-dimensional parameters of the parabolic velocity 
boundary condition problem. 


Parameter Value 

a 1 and 10 

/i 1 and 0.001 

P 1 

b(x) [1, 1] 

L 1 


Table 4. Pressure slab problem: Non-dimensional parameters used in the problem. 


Parameter Value 


a 

1 

P 

0.001 

P 

1 

Pinj 

5 and 7.5 

Patm 

1 

L 

1 

W 

0.2 


bottom sides of the boundary. The traction is prescribed on the left side (i.e., tP(x) = —pinjn(x) 
on r^) and on the right side (i.e., tP(x) = —patmn(x) on P^). We shall use this test problem to 
assess the accuracy of numerical solutions with respect to the reciprocal relation. Figure 12 depicts 
the variation of the error in the reciprocal relation with mesh refinement for the triangular and 
quadrilateral grids. The error in the reciprocal relations under Darcy equations is very close to 
zero for all the meshes. The error under Darcy-Brinkman equations decrease uniformly with mesh 
refinement. All the numerical results are in accordance with the theoretical predictions. 

4.5. Pressure driven problem. The domain is a square with L = 1. A pictorial descrip¬ 
tion of the problem is given in Figure 13. The traction is prescribed on the left boundary (i.e., 
tP(x) = —pinjn(x) on F^) and on the middle of the right boundary (i.e., tP(x) = —patmn(x) on 
F 2 ). Elsewhere, homogeneous velocity boundary condition is enforced. Due to the prescription of 
the traction on a part of the boundary, this problem is incompatible with the minimum dissipa¬ 
tion theorem but is compatible with the reciprocal relation and total mechanical power theorems. 
Herein, we present the results for the reciprocal relation. The non-dimensional parameters used in 
the numerical simulation are provided in Table 5. Figure 14 shows the variation of the error in the 
reciprocal relation with mesh refinement for triangular and quadrilateral meshes. The convergence 
under uniform structured meshes is uniformly to zero which is expected by theorem. 

4.6. Vorticity results. The maximum principle given by Theorem 4 can be used to assess 
the accuracy of numerical solutions to the Darcy-Brinkman equations by checking whether the non¬ 
negative maximum and non-positive minimum of the vorticity occur on the boundary. Figure 15 
shows that the maximum principle is satisfied for various problems under the steady-state Darcy- 
Brinkman equations. 
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Table 5. Pressure driven problem: Non-dimensional parameters used in the problem. 


Parameter 

Value 

a 

1 

tJ- 

1 and 0.001 

P 

1 

Pin] 

5 and 7.5 

Patra 

1 

L 

1 


5. SYNTHETIC RESERVOIR DATA: MARMOUSI DATASET 

We will now solve an idealized reservoir problem using a popular synthetic dataset - the so- 
called (smooth) Marmousi dataset [Versteeg and Grau, 1990; Versteeg and Lailly, 1991; Versteeg, 
1993; Klimes, 2014; Benamou, 2014], The dataset provides spatially varying speed of sound on a 
384 X 122 grid. We have assumed that the permeability scales linearly with the values provided by 
the dataset. This is just an arbitrary choice to generate a heterogeneous dataset for permeability. 
However, it should be noted that the conclusions that will be drawn here will be valid even if one 
uses another dataset for the permeability. Figure 16 shows the contours of Marmousi dataset. 

The boundary value problem of the reservoir is pictorially described in Figure 17. This compu¬ 
tational domain has been employed in some recent works (e.g., [Nakshatrala and Rajagopal, 2011]). 
All these works have assumed homogeneous medium properties, and did not use a reservoir data like 
the Marmousi dataset. Moreover, these studies did not address the use of a posteriori techniques 
to access numerical accuracy which is the main focus of the current paper. The parameters used in 
this problem are provided in Table 6. Figure 18 shows that the errors in the reciprocal relation are 
larger under uniform structured meshes than under adaptive meshes. This is due to the fact that 
uniform structured meshes suffer from pollution errors due to the singularity near the production 
well. The variation of the error in the reciprocal relation with mesh refinement provides guidelines 
on how much rehnement is required to obtain solution of some desired accuracy especially in those 
cases where there are no analytical solutions. Figure 19 shows the magnitude of the vorticity, and 
one can see from the figure that the maximum principle for the vorticity is satisfied. 

Table 6. Synthetic reservoir problem: Non-dimensional parameters used in the problem. 


Parameter 

Value 

P 

0.001 

P 

1 

Pin] 

5 and 7.5 

Patm 

1 

L 

384 

H 

384/2 

W 

384x0.1 

k 

Marmousi dataset/L^ 
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6. TRANSIENT CASE 


6.1. Governing equations. Let us denote the time interval of interest by X, and the time by 
t £l. The unsteady governing equations under the Darcy model take the following form: 


d'v 

p— -h av -h gradjp] = pb(x, t) 

in D X X 

(6.1a) 

div[v] = 0 

in D X X 

(6.1b) 

v(x,t) • n(x) = i;n(x,t) 

on X X 

(6.1c) 

p(x,t) =po(x,t) 

on r* X X 

(6.1d) 

v(x,t = 0) = vo(x) 

in D 

(6.1e) 


Of course, the convective term grad[v]v is neglected. We assume that the coefficient of viscosity of 
the fluid and the permeability of the porous solid to be constants, and hence the drag coefficient is 
constant. We further assume that the density is homogeneous, and the body force is assumed to 
be conservative. The unsteady Darcy-Brinkman equations can be written as follows: 


p— + av + grad[p] — div[2^D] = /3b(x, t) in D x X (6.2a) 

div[v] =0 in D X X (6.2b) 

v(x, t) = vP(x, t) on r’'X X (6.2c) 

Tn(x) = tP(x, t) on r* X X (6.2d) 

v(x, t = 0) = vo(x) in D (6.2e) 


where vP(x, t) is the prescribed velocity vector, and tP(x, t) is the prescribed traction. 


6.2. Mathematical properties. Under unsteady Darcy equations, the vorticity satisfies the 
following equation: 

du) 

p— + au = 0 (6.3) 

which is an ordinary differential equation at each spatial point. The solution takes the following 
form: 


cj(x, t) = a;o(x) exp 


at 

P 


(6.4) 


This means that the decay of vorticity should be exponential. One can check whether the numerical 
solutions exhibit this trend by plotting in log scale the each component of the vorticity with respect 
to time should, and this should be a straight line with slope —ajp. If this trend is not satisfied 
for a given mesh and given time-step, does refining the grid spacing and the time-step improve the 
trend? Of course, one can check whether the vorticity goes to zero for large times. The vorticity 
under unsteady Darcy-Brinkman equations satisfies the following equation: 

du) 

P~^ + “ /xAo; = 0 in D X X (6.5) 

where A is the Laplacian operator. The above equation is a homogeneous linear parabolic partial 
differential equation, which is known to satisfy a maximum principle [Pao, 1993]. This implies that 
both the maximum and the minimum will occur either in the initial condition or on the boundary. 
One can check whether numerical solutions satisfy the aforementioned maximum principle. Also, 
whether refining the grid spacing and time-steps affect the performance of numerical solutions with 
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Table 7. Pressure slab problem: Initial condition used in the problem. 


Parameter Value 

P Patm 

Vx = Vy sin( 7 rx/IT) sin( 7 ry/L) 


respect to this metric. One can devise any test problem as long as the following assumptions are 
met: (i) p is constant, (ii) permeability is homogeneous, and (iii) the body force is conservative. 

6.3. Representative numerical results. The theoretical results presented in this section 
are corroborated numerically in Figures 20-21. The vorticity results for the pressure slab problem 
under the transient Darcy model are provided in Figures 20. (Recall that Figure 11 provided 
a pictorial description of the pressure slab problem.) The non-dimensional parameters used in 
the numerical simulation are presented in Table 4, and the initial condition is provided in Table 
7. Figure 20 shows that log(^) is a straight line with slope equal to ^ under the mesh and time 
refinements for various spatial points in the solution of the transient Darcy equations (for the current 
problem ^ = —1). Figure 21 verifies the maximum principle for vorticity for various test problems 
under transient Darcy-Brinkman equations. In all the cases, the non-negative maximum and non¬ 
positive minimum of the vorticity occur on the boundary, which agree with the mathematical theory 
presented above. 


7. CONCLUDING REMARKS 

We presented a new methodology for solution verification which is independent of utilized 
numerical formulation. We called it mechanics-based solution verification. To employ this approach 
we developed four important properties (minimum dissipation theorem, minimum total mechanical 
power, reciprocal relation, and maximum principle for vorticity) that the solutions under the Darcy 
and the Darcy-Brinkman models satisfy. However, the proposed technique is not merely restricted 
to porous media type problems and one can easily extend it to other models. We also presented 
various test problems can serve as benchmark problems to verify numerical implementation of 
solvers for the Darcy and Darcy-Brinkman models. Results showed that the proposed technique 
can be effectively used to assess the accuracy of numerical solutions, identify numerical pollution, 
and check the performance of adaptive mesh. An attractive feature is that these properties can 
be verified for any given problem (i.e., needed not be one of the benchmark problems), and for 
any computational domain. For example, if the problem involves prescribing velocity boundary 
condition on the entire boundary, then one plots the dissipation with respect to mesh refinement. 
Some of the main conclusions are: 

(a) If the numerical formulation is converging, the dissipation, total mechanical power, and recipro¬ 
cal relation should decrease with mesh refinement. If this does not occur, one needs to suspect 
that there are singularities in the solutions or that the numerical formulation does not perform 
well with respect to the local mass balance property. 

(b) It has been shown that the minimum dissipation, minimum total mechanical power, and re¬ 
ciprocal relation theorems can be utilized to identify pollution errors in numerical solutions. 
The theorems can also be used to assess whether a given type of mesh will be able to resolve 
singularities in the solution. This can be assessed by creating a series of hierarchical meshes and 
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plotting for example the dissipation, with respect to characteristic mesh size. A given type of 
mesh will resolve singularities in the solution and will not be affected by pollution errors if the 
total dissipation decreases uniformly and reaches a plateau with a hierarchical mesh refinement. 

(c) The non-negative maximum vorticity and the non-positive minimum vorticity under Darcy- 
Brinkman equations with homogeneous isotropic medium properties should occur on the bound¬ 
ary. 

(d) Log plot of vorticity with respect to time, under transient Darcy models, is a line with slope 

_ a 

P ' 

The proposed a posteriori techniques can be invaluable additions to the usual repertoire of methods 
for verification - the method of exact solutions (MES) and the method of manufactured solutions 
(MMS). 


APPENDIX: UNIQUENESS OF SOLUTION 

For completeness and pedagogical value, we now show that the uniqueness of solution under 
Darcy-Brinkman equations is a direct consequence of the minimum dissipation inequality. One can 
construct a proof for the uniqueness of solution under Darcy equations on similar lines. 

Theorem 5 (Uniqueness theorem). The solution to Darcy-Brinkman equations (2.4a)-(2.4d) is 
unique up to an arbitrary constant for the pressure given b(x), vP(x) and tP(x). 

Proof. On the contrary assume that {vi(x),pi(x)} and {v 2 (x),p 2 (x)} are two solutions to 
Darcy-Brinkman equations for the prescribed data. Let us consider the following quantity: 

X := / a(x) (vi(x) - V 2 (x)) • (vi(x) - V 2 (x)) dQ-\- [ 2//(Di(x) - D 2 (x)) • (Di(x) - D 2 (x)) dD 
J fi Jq 

Noting that div[vi] = 0 and div[v 2 ] = 0, the second integral can be simplified as follows: 

[ 2^(Di(x)-D2(x))-(Di(x)-D2(x)) dQ= [ (Di(x) - D 2 (x)) • (Ti(x) - T 2 (x)) dD 
Jq J fi 

= f grad [vi(x) - V 2 (x)] • (Ti(x) - T 2 (x)) dD 
Jn 

In obtaining the above equation, we have used the fact that Ti(x) and T 2 (x) are symmetric. Using 
Green’s identity, the above equation can be written as follows: 

[ 2|U (Di(x) - D 2 (x)) • (Di(x) - D 2 (x)) dD = [ (vi(x) - V 2 (x)) • (Ti(x) - T 2 (x))n(x) dT 
Jn Jan 

- [ (vi(x) - V 2 (x)) • (div[Ti] - div[T 2 ]) dD 

Jn 

Using boundary conditions and the balance of linear momentum, we get the following: 

[ 2p, (Di(x) - D 2 (x)) • (Di(x) - D 2 (x)) dD = - / a(x) (vi(x) - V 2 (x)) • (vi(x) - V 2 (x)) dD 
Jn J fi 

This implies that X = 0. Since a(x) > 0 Vx G D and /U > 0 in D, one can conclude that 

vi(x)=V 2 (x) Vx G D 
Di(x) = D 2 (x) Vx G D 
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(7.1a) 

(7.1b) 



That is, the velocity and symmetric part of velocity vector field are unique. Using the equation for 
the balance of linear momentum (2.4a) and the fact that the velocity vector field is unique, one can 
obtain the following equation: 

grad [pi(x) — P 2 (x)] = 0 Vx G 14 (7.2) 

This further implies that 

pi(x)-p 2 (x) =po Vx G n (7.3) 

where po is an arbitrary constant. This completes the proof. □ 
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O Location of pressure unknown • Location of velocity unknown 

(c) Q2P1 interpolation for nine-node quadrilat- (d) P2P1 interpolation for six-node triangular 

eral (Q9) element (T6) element 


Figure 1. This figure shows the typical structured finite element meshes employed 
in the current study. We use Q2P1 and P2P1 mixed interpolations to approximate 
the unknowns (i.e., second-order interpolation for the velocity field, and first-order 
interpolation for the pressure field). 
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vP(x) = 0 or v„(x) = 0 


Figure 2. Body force problem: The computational domain is a square with L = 1. 
The prescribed conservative body force is /9b(x) = 10 x [sin(7rx), cos(7r?/)]. Homoge¬ 
neous velocity is enforced on the entire boundary (i.e., F^ = dfl). That is, vP(x) = 0 
for Darcy-Brinkman equations and u„(x) = 0 for Darcy equations. 
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(a) Darcy model with ^ = 1. 


(b) Darcy-Brinkman model with /r = 0.001. 


Figure 3. Body force problem: This figure shows the variation of the dissipation 
and the total mechanical power with mesh refinement under the Darcy and Darcy- 
Brinkman models using quadrilateral and triangular elements. The parameters used 
in this problem are provided in Table 1. The figure show that the total mechanical 
power and the dissipation decrease uniformly with mesh refinement, which are in 
accordance with Theorems 1 and 2. 




(b) Darcy-Brinkman model with = 0.001. 


Figure 4. Body force problem: This figure shows the variation of Ereciprocai with 
mesh refinement for the Darcy and Darcy-Brinkman models using quadrilateral and 
triangular elements. The parameters used in this problem is provided in Table 1. 
The figure show that the numerical error in the reciprocal relation decrease uniformly 
to zero with mesh refinement for this test problem. 
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Figure 5. Lid-driven cavity problem: The computational domain is a unit square. 
Velocity boundary condition is prescribed on the entire boundary (i.e., F^ = dQ). 
The prescribed velocity on the top side is Ua; = 1 and Vy = 0, and the prescribed 
velocity on the remaining sides of the boundary is zero (i.e., vP(x) = 0). Note that 
this test problem is not compatible with Darcy equations. 


Results for Lid-driven cavity problem 
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(a) Uniform structured mesh with 
Ijh = 12 . 



(c) Dissipation vs. h under a hierarchy of uni¬ 
form structured meshes. 
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(b) Adaptive mesh near the top left corner 
with l//i = 12 elsewhere. 



(d) Dissipation vs. h under a hierarchy of 
adaptive meshes. 



Figure 6. Lid-driven cavity problem: The top figures show a uniform structured 
mesh, and an adaptive mesh near the top corners. The bottom figures show the 
variation of total dissipation with mesh refinement for quadrilateral and triangular 
finite elements. The parameters in this problem are provided in Table 2. One should 
note that the solution exhibits singularities at the top corners. It is evident that 
the convergence is uniform under both structured and adaptive meshes. However, 
the dissipation under the structured mesh increased with mesh refinements, whereas 
the trend is reversed under the adaptive mesh. This is because the adaptive mesh 
resolves the singularities and removes the associated pollution error. This figure 
corroborates one of the main conclusions of the paper that the minimum dissipation 
theorem can he utilized to identify pollution errors in numerical solutions. 
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v(x) = 0 



Figure 7. Pipe bend problem (velocity boundary condition): A pictorial description of 
the problem. The computational domain O is a unit square. The velocity boundary 
condition is prescribed on the entire boundary (i.e., T’^ = d^l). As indicated in the 
figure, a parabolic velocity profile with = 1 is prescribed on a segment of the 

left side of the boundary. Similarly, a parabolic velocity profile with = 1 is 

prescribed on a segment of the bottom side of the boundary. Homogeneous velocity 
is enforced on the remaining parts of the boundary. 
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(b) Total mechanical power for a = 1. 


Figure 8. Pipe bend problem (velocity boundary condition): The figure shows that 
the dissipation and total mechanical power decrease uniformly with mesh refinement 
for the Darcy and Darcy-Brinkman models. The results are presented for both 
quadrilateral and triangular elements. The parameters used in this problem are 
provided in Table 3 (;U = 1 for Darcy and /U = 0.001 for Darcy-Brinkman models). 
Note that the total mechanical power is defined in equation (3.3). 
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v(x) = 0 



Figure 9. Pipe bend problem (velocity-pressure boundary condition): This figure 
presents a pictorial description of the problem. The computational domain Q is 
a square with L = 1. The traction boundary condition is tP(x) = — patmn(x) on T*. 
A parabolic velocity profile with = 1 is prescribed on a segment of the left 

side of the boundary, as indicated in the figure. Elsewhere, the velocity is assumed 
to be zero (i.e., vP(x) = 0 for Darcy-Brinkman equations and Un(x) = 0 for Darcy 
equations). Note that the corners at the outlet are re-entrant corners. 
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(a) Uniform structured mesh with l/h = 
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(b) Adaptive mesh near the outlet with 
l/h = 15 elsewhere. 
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(c) Darcy-Brinkman model under the uniform (d) Darcy-Brinkman model under the adaptive 
structured mesh, fi = 0.001. mesh, jj, — 0.001. 


Figure 10. Pipe bend problem (velocity-pressure boundary condition); The top fig¬ 
ures show the uniform structured mesh and the adaptive mesh near the corners of 
the outlet (i.e., F*). The bottom figures show the variation of the total mechanical 
power with mesh refinement under quadrilateral and triangular elements. The pa¬ 
rameters used in this problem are provided in Table 3. The total mechanical power 
increased uniformly with mesh refinement under the structured mesh, while the total 
mechanical power decreased uniformly and reached a plateau with mesh refinement 
under the adaptive mesh. The numerical results clearly demonstrate that the total 
mechanical theorem can be utilized to identify pollution errors due to singularities 
and to assess whether a particular type of computational mesh is suitable for a given 
problem. 
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Figure 11. Pressure slab problem: A pictorial description of the problem. The 
computational domain is aW x L rectangle. The traction is prescribed on the left 
side of the boundary (i.e., tP(x) = — pinjn(x) on T^) and on the right side (i.e., 
tP(x) = —patmn(x) on r^). Elsewhere, homogeneous velocity is enforced. 



Figure 12. Pressure slab problem: The figure shows the variation of eredprocai with 
mesh refinement for Darcy and Darcy-Brinkman equations using quadrilateral and 
triangular grids. The parameters in this problem are provided in Table 4. One 
can see the error in the numerical results with respect to the reciprocal relation for 
Darcy equations is very close to zero for all meshes. The corresponding error under 
Darcy-Brinkman equations decreases uniformly with mesh refinement. 
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Figure 13. Pressure-driven problem: A pictorial description of the problem. The 
computational domain is a square with L = 1. The traction is prescribed on the 
left side of the boundary (i.e., tP(x) = —pinjn(x) on T^) and on the middle of right 
side (i.e., tP(x) = —patmn(x) on T^). Elsewhere, homogeneous velocity is enforced. 



Figure 14. Pressure-driven problem: Figure shows the variation of ^reciprocal under 
the Darcy-Brinkman models using quadrilateral and triangular finite elements. The 
parameters used in this problem are provided in Table 5. Clearly, the reciprocal 
theorem can be utilized to obtain information on the numerical performance of the 
problem that does not possess an analytical solution. 
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(a) Lid-driven cavity, regular mesh. LUmax 
35.641 and uJmin = 0. 


(b) Lid-driven cavity, adaptive mesh. LUmax 
778 and = 0. 
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(c) Body force problem, ujrnax = 0.1923 and oJmin. = (d) Slab problem. oJmax = 561.52 and 

4.8073 X 10"®. ujn^in = 3.1957 x 10"®. 


Figure 15. The figure verifies the maximum principle for the vorticity for var¬ 
ious two-dimensional problems under the Darcy-Brinkman model. We employed 
quadrilateral elements with size 1/h = 20. The numerical results corroborate the 
theoretical predictions given in Theorem 4 for all the test problems we have consid¬ 
ered. 
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Figure 16. Synthetic reservoir problem: This figure shows the contours of (smooth) 
Marmousi dataset [Benamou, 2014], The dataset provides values on a 384 x 122 
grid which we have scaled to our rectangular computational domain oi L = 2 and 
H = 1. This spatially heterogeneous dataset is widely used as a benchmark dataset 
in reservoir modeling. 


Results for Synthetic reservoir problem 
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Figure 17. Synthetic reservoir problem: This figure provides a pictorial description 
of the test problem. The domain of the problem is a rectangle of size H x L. 
The injection pressure is prescribed on the left and right boundaries (i.e., tP(x) = 
—Pinjn(x) on r^), and the atmosphere pressure is prescribed on the middle of top 
side (i.e., tP(x) = —patmS(x) on T^). The prescribed velocity on the remaining 
parts of the boundary is zero (i.e., vP(x) = 0 for the Darcy-Brinkman model and 
Un(x) = 0 for the Darcy model). 
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(c) Variation of Sreciprocai under a hierarchy of uni- (d) Variation of Sreciprocai under a hierarchy of adaptive 
form structured meshes. meshes. 


Figure 18. Synthetic reservoir problem: The figures show the variation of ^reciprocal 
with h for the Darcy and Darcy-Brinkman models using structured and adaptive 
meshes. The parameters in this problem are provided in Table 6. The figures 
clearly show that the errors in the reciprocal relation are larger under the uniform 
structured meshes. This can be attributed to pollution errors due to singularities 
at the corners of the production wells. 
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Figure 19. Synthetic reservoir problem: The figure shows the magnitude of vor- 
ticity for two-dimensional Darcy-Brinkman equations using quadrilateral elements. 
The parameters used in this problem are provided in Table 6. This figure verifies 
the maximum principle for the vorticity for the synthetic reservoir problem using 
Marmousi dataset. 
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(b) Initial vorticity for quadrilateral ele¬ 
ment size of l//i = 20. 




(c) Mesh refinement for quadrilateral elements, at (d) Time refinement for 1/h = 20 at x = 

X = (0.05,0.95) and dt = 0.01. (0.15,0.35). 


Figure 20. Pressure slab problem: The figure verifies the theoretical results for 
the vorticity under transient Darcy equations. The results show that the slope of 
log(^) for various spatial points in the computational domain are close to — ^ = 
— 1. The hgure also shows the slope is follows the theoretical prediction for various 
mesh rehnements h and time-steps dt. The non-dimensional parameters used in this 
numerical experiment are provided in Tables 4 and 7. 
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(a) Lid-driven cavity, structured mesh. 


(b) Lid-driven cavity, adaptive mesh. 
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(c) Body force problem. 


(d) Pressure slab problem. 


Figure 21. The figure verifies the maximum principle for the vorticity under tran¬ 
sient Darcy-Brinkman equations for various two-dimensional problems. The numer¬ 
ical results satisfy the maximum principle for all the test problems. The results are 
reported at t = 1 and for time-step dt = 0.5 (which is chosen arbitrarily). The 
results are presented for the adaptive and structured meshes based on quadrilat¬ 
eral elements with size 1/h = 10. We employed the backward Euler time-stepping 
scheme in the numerical experiment. 
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